#from gtk._gtk import Orientable
from tvtk.plugins.scene.ui.actions import YMinusView
PREFIX_U = 'u_r_'
PREFIX_V = 'v_r_'
PREFIX_RHO = 'rho_r_'
PATH = '/home/behollis/DATA/in/ts00050/'
DPATH = '/home/behollis/DATA/out/lockExSt/ts00050/'
#ODPATH = '/home/behollis/Dropbox/eurovis2015shortpaperCode/oslines/'
OUTPUT_DATA_DIR = '/home/behollis/Dropbox/eurovis2015shortpaperCode/oslines/'

import math
import sys
import loadStreamlines
import h5py
from loadStreamlines import readStreamlines
import interpProbDensity
import matplotlib.pyplot as plt
import numpy as np
import matplotlib
from skimage.draw import line_aa

prob = None
slines = None


def openH5pyFiles():
     #print 'reading hdf5 file...'
    #backward integrated...
    '''
    f0 = h5py.File(DPATH+'0.0.hdf5', 'r')
    f1 = h5py.File(DPATH+'0.1.hdf5', 'r')
    f2 = h5py.File(DPATH+'0.2.hdf5', 'r')
    f3 = h5py.File(DPATH+'0.3.hdf5', 'r')
    f4 = h5py.File(DPATH+'0.4.hdf5', 'r')
    f5 = h5py.File(DPATH+'0.5.hdf5', 'r')
    '''
    
    f0 = h5py.File(DPATH+'1.6.hdf5', 'r')
    f1 = h5py.File(DPATH+'1.7.hdf5', 'r')
    f2 = h5py.File(DPATH+'1.8.hdf5', 'r')
    f3 = h5py.File(DPATH+'1.9.hdf5', 'r')
    f4 = h5py.File(DPATH+'1.10.hdf5', 'r')
    f5 = h5py.File(DPATH+'1.11.hdf5', 'r')
    
    #f6 = h5py.File(DPATH+'1.6.hdf5', 'r')
    #f7 = h5py.File(DPATH+'1.7.hdf5', 'r')
    #f8 = h5py.File(DPATH+'1.8.hdf5', 'r')
    #f9 = h5py.File(DPATH+'1.9.hdf5', 'r')
    #f10 = h5py.File(DPATH+'1.10.hdf5', 'r')
    #f11 = h5py.File(DPATH+'1.11.hdf5', 'r')
    
    #f12 = h5py.File(DPATH+'2.12.hdf5', 'r')
    #f13 = h5py.File(DPATH+'2.13.hdf5', 'r')
    #f14 = h5py.File(DPATH+'2.14.hdf5', 'r')
    #f15 = h5py.File(DPATH+'2.15.hdf5', 'r')
    #f16 = h5py.File(DPATH+'2.16.hdf5', 'r')
    #f17 = h5py.File(DPATH+'2.17.hdf5', 'r')
    
    '''
    f0f = h5py.File(DPATH+'0.0.forward.hdf5', 'r')
    #f1f = h5py.File(DPATH+'0.1.forward.hdf5', 'r')
    f2f = h5py.File(DPATH+'0.2.forward.hdf5', 'r')
    #f3f = h5py.File(DPATH+'0.3.forward.hdf5', 'r')
    f4f = h5py.File(DPATH+'0.4.forward.hdf5', 'r')
    #f5f = h5py.File(DPATH+'0.5.forward.hdf5', 'r')
    
    f6f = h5py.File(DPATH+'1.6.forward.hdf5', 'r')
    #f7f = h5py.File(DPATH+'1.7.forward.hdf5', 'r')
    f8f = h5py.File(DPATH+'1.8.forward.hdf5', 'r')
    #f9f = h5py.File(DPATH+'1.9.forward.hdf5', 'r')
    f10f = h5py.File(DPATH+'1.10.forward.hdf5', 'r')
    #f11f = h5py.File(DPATH+'1.11.forward.hdf5', 'r')
    
    f12f = h5py.File(DPATH+'2.12.forward.hdf5', 'r')
    #f13f = h5py.File(DPATH+'2.13.forward.hdf5', 'r')
    f14f = h5py.File(DPATH+'2.14.forward.hdf5', 'r')
    #f15f = h5py.File(DPATH+'2.15.forward.hdf5', 'r')
    f16f = h5py.File(DPATH+'2.16.forward.hdf5', 'r')
    #f17f = h5py.File(DPATH+'2.17.forward.hdf5', 'r')
    '''
    
    print 'finished reading h5py file!'
    
    #iles = ( f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11, f12, f13, f14, f15, f16, f17 )
    #files_f = ( f0f,f1f,f2f,f3f,f4f,f5f,f6f,f7f,f8f,f9f,f10f,f11f, f12f, f13f, f14f, f15f, f16f, f17f )
    
    files = ( f0, f1, f2, f3, f4, f5 )
    #files_f = ( f0f,f1f,f2f,f3f,f4f,f5f,f6f,f7f,f8f,f9f,f10f,f11f, f12f, f13f, f14f, f15f, f16f, f17f )
         
    return files#, files_f

def closeH5pyFiles(files):
    for f in files:
        f.close()
        
def readStreamlines(x, y, b_files, f_files, dinteg = 'backward', max_per_file=100):
    ''' loads streamlines thru seed with spatial extents / member extents '''
    
    slines = dict()
    xlst =list()
    ylst = list()
    zlst = list()
    
    cnt = 0
   
    for f in b_files:
        #print 'loading file'+str(file)
        
        groups = f.keys() #use only f for keys, but they should be the same as for ff
        for mem in groups:
            
            if cnt > max_per_file:
                cnt = 0
                break
            
            dir = str(mem).zfill(4) + '/x' + str(x).zfill(3) + '/y' + str(y).zfill(3)
        
            try:
                if dinteg == 'backward':
                    xlst = list(f[dir][0])
                    ylst = list(f[dir][1])
                    zlst = 0.0 * len(list(f[dir][0]))
                '''
                elif dinteg == 'forward':
                    xlst = list(ff[dir][0])
                    ylst = list(ff[dir][1])
                    zlst = 0.0 * len(list(ff[dir][0]))
                else: # 'both'
                    xlst_b = list(f[dir][0][1:])
                    xlst_b.reverse()
                    
                    ylst_b = list(f[dir][1][1:])
                    ylst_b.reverse()
                    
                    zlst_b = [0.0] * len( xlst_b )
                    #zlst_b.reverse()
                    
                    xlst_f = list(ff[dir][0])
                    ylst_f = list(ff[dir][1])
                    zlst_f = [0.0] * len( xlst_f )
                    
                    xlst = xlst_b + xlst_f
                    ylst = ylst_b + ylst_f
                    zlst = zlst_b + zlst_f
                '''
                    
                slines[mem] = (xlst, ylst, zlst)
                
                cnt += 1
                    
            except:
                print 'trying to read MEMBER or SEED that is not in this hdf5 file...'
                
    return slines

def calcStreamlineProbs2(seed_pt, sl_dict, MEM):
    
    slines_prob = dict()
    
    # gather all of the vector fields
    #print len(sl_dict.keys())
    
    for midx in sl_dict.keys():
        #print midx
        vfield = interpProbDensity.readMember( int(midx[3:]) )
        interpProbDensity.gUv[midx] = vfield
    
    probs_list = list()#contains all prob densities
    
#    for slidx, sl in enumerate(slines):
    for slkey in sl_dict.keys():
        sl = sl_dict[slkey]
        #tot = 0
        slines_prob[slkey] = [ 0.0, [] ]
        for ptidx in range(0, len(sl[0])):
            
            #print 'streamline length: ' + str( len(sl[0]) )
            #print 'for sl key: ' + str(slkey)
            
            pt = [ sl[0][ptidx], sl[1][ptidx] ]
            
            kde, min_u, max_u, min_v, max_v = interpProbDensity.interpVelFromEnsemble(pt[0:2], MEM)
            #plotDistro( kde, (min_u,max_u), (min_v,max_v), '') 
            #band_width = MEM**(-1./(DIM+4))
            EPS_FACTOR = 0.1
            epsx = EPS_FACTOR * math.fabs(max_u - min_u)
            epsy = EPS_FACTOR * math.fabs(max_v - min_v)
            
            u = interpProbDensity.gUv[slkey][0][pt[0],pt[1]]
            v = interpProbDensity.gUv[slkey][1][pt[0],pt[1]]
            low_bounds = [u - epsx, v - epsy]
            high_bounds = [u + epsx, v + epsy]
            prob = kde.integrate_box( low_bounds, high_bounds )
            
            if prob < 0.0:
                prob = 0
            
            #SCALE = 100
            #prob *= SCALE
            
           
            slines_prob[slkey][0] += prob  
            
            #print prob
            #print slines_prob[slkey][0]
            
            #tot += prob #Possibly make this a Markov chain (however densities are very low)
                
            #print 'for pt: ' + str(ptidx)
            #print 'prob' + str(prob)
            #print 'tot' + str(tot)
                
            slines_prob[slkey][1].append(float(prob))
            
            
            probs_list.append(prob)
            
                
            #print('streamline {0}: {1}'.format(slkey, tot))
            
        
    return slines_prob, probs_list

from operator import itemgetter
import random

def renderStandardPlot(top=None):
    ''' 
    Rank each streamline based on total KDE prob density #
    Render all for comparison
    '''
    
    fig = plt.figure(frameon=False)
    plt.box('off')
    plt.imshow(im)
    
    plt.tick_params(\
    axis='both',          # changes apply to the x-axis
    which='both',      # both major and minor ticks are affected
    bottom='off',
    right='off',
    left='off',      # ticks along the bottom edge are off
    top='off',         # ticks along the top edge are off
    labelbottom='off', labelleft='off') # labels along the bottom edge are off
    
    '''
    ax = fig.gca(projection='3d')
    axes = fig.add_subplot(1, 1, 1, axisbg='white')
    axes.set_aspect('equal', 'datalim')
    '''
    # render streamlines
    for slkey in sl_dict.keys():
        sl = sl_dict[slkey]
        if top is not None and slkey in top:
            sl = np.asarray(sl)
            flip = np.zeros( len(sl[1]) ); flip.fill( im.shape[1] )
            plt.plot( IMGSCALEX*np.asarray( sl[0] ), \
                     flip - IMGSCALEY*np.asarray(sl[1]),\
                     color=(0.8, 0.0, 0.0), linewidth = 1.0 )
     
    # draw seed loc  
    plt.plot([ IMGSCALEX*sx ], [ im.shape[1] - IMGSCALEY*sy ], '+', color='black',linewidth = LW, ms=MS )  
    plt.xlim([xmin, xmax])
    plt.ylim([ymin, ymax])
    #plt.show()
    #lt.gcf().set_size_inches(,size/float(dpi)))
    plt.savefig('../img/all.png')

def renderOverallLikelihoods(top=None):
    from numpy import linalg
    
    fig = plt.figure()
    plt.box('off')
    plt.imshow(im)
    
    plt.tick_params(\
    axis='both',          # changes apply to the x-axis
    which='both',      # both major and minor ticks are affected
    bottom='off',
    right='off',
    left='off',      # ticks along the bottom edge are off
    top='off',         # ticks along the top edge are off
    labelbottom='off', labelleft='off') # labels along the bottom edge are off
    
    # Render with likelihoods
    for tup in toplist:
        slkey = tup[0]
        if top is not None and slkey in top:
            sl = sl_dict[slkey]
            tot = tup[1]
            clr = mapper2.to_rgba(x=tot) 
            
            alpha = tot / max2 #normalize magnitude for alpha                             
            flip = np.zeros( len(sl[1]) ); flip.fill( im.shape[1] )
            
            plt.plot( IMGSCALEX*np.asarray( sl[0] ), \
                     flip - IMGSCALEY*np.asarray(sl[1]),\
                     color=(clr[0], clr[1], clr[2], alpha), linewidth = 1.0 )
            
    plt.plot([ IMGSCALEX*sx ], [ im.shape[1] - IMGSCALEY*sy ], '+', color='black',linewidth = LW, ms=MS )
    plt.xlim([xmin, xmax])
    plt.ylim([ymin, ymax])  
    #plt.colorbar(mapper2)   
    #plt.show() 
    plt.savefig('../img/topProb.png')    

def renderIndividualSlines():
    # Rendering streamlines by velocity prob density per integration step
    
    '''
    http://matplotlib.org/examples/pylab_examples/multicolored_line.html
    http://stackoverflow.com/questions/19390895/matplotlib-plot-with-variable-line-width
    '''
    
    from matplotlib.collections import LineCollection
    
    for slkey in sl_dict.keys():
        
        fig = plt.figure()
        #ax = fig.gca(projection='3d')
        axes = fig.add_subplot(1, 1, 1, axisbg='white')
        axes.set_aspect('equal', 'datalim')
        # render streamlines
        
        sl = sl_dict[slkey]
       
        try:
            l = list(slines_prob[slkey][1])
            #print 'length of list l: ' + str( len(l) )
            
            min_= min(l)
            max_ = max(l)
        except:
            #print 'streamline not more than one integration step'
            continue
        
        # draw each segment a color depending on prob density
        steps = len( sl[0] ) - 1
        
        x = list()
        y = list()
        lws = list()
        clrs = list()
        
        for idx in range(0, steps):
           
            #normalize by shifting / scaling to [0,1]
            scale = (slines_prob[slkey][1][idx] - min_0) / math.fabs( max_0 - min_0) 
        
                
            prob = slines_prob[slkey][1][idx]
            print prob
            '''
            if scale < 0:
                scale = 0.0
            print scale
            '''
            clr = mapper.to_rgba(x=prob) 
            print clr  
            #clr = (1.0*scale, 1.0*scale, 1.0*scale, 1.0)
            x0 = sl[0][idx]; x1 = sl[0][idx+1]
            y0 = sl[1][idx]; y1 = sl[1][idx+1]
            SCALE_W = 3.0
            MIN_W = 1.0
            #axes.plot([x0, x1 ], [y0, y1], color=clr, linewidth = scale*SCALE_W + MIN_W )
            x.append(x0); x.append(x1)
            y.append(y0), y.append(y1)
            lws.append(scale*SCALE_W + MIN_W)
            clrs.append(clr)
            
        points = np.array([x,y]).T.reshape(-1,1,2)
        segments = np.concatenate([points[:-1], points[1:]], axis=1)
        lc = LineCollection(segments, linewidths=lws, colors=clrs)
        axes.add_collection(lc)
        axes.plot([ sx ], [ sy ], '+', color='black', linewidth=LW , ms=MS) # draw seed loc
            
        #plt.colorbar(mapper)
        
        #plt.xlim(x0=50)
        #plt.ylim(x0=50)
        #plt.xlim([50,100])
        #plt.ylim([50,100])
        plt.savefig('../img/sline{0}.png'.format(slkey))

def renderColorbars():
    #get images of colorbars
    # make a color map of fixed colors
    fig, ax = plt.subplots()
    fig.colorbar(mapper, orientation='horizontal', ticks=[min_0, max_0] )
    plt.savefig('../img/colorbar1.png')
    fig, ax = plt.subplots()
    fig.colorbar(mapper2, orientation='horizontal', ticks=[min2, max2])
    plt.savefig('../img/colorbar2.png')
    
def get_line(x1, y1, x2, y2):
    '''
    http://www.roguebasin.com/index.php?title=Bresenham%27s_Line_Algorithm
    http://www.cs.helsinki.fi/group/goa/mallinnus/lines/bresenh.html
    '''
    
    points = []
    issteep = abs(y2-y1) > abs(x2-x1)
    if issteep:
        x1, y1 = y1, x1
        x2, y2 = y2, x2
    rev = False
    if x1 > x2:
        x1, x2 = x2, x1
        y1, y2 = y2, y1
        rev = True
    deltax = x2 - x1
    deltay = abs(y2-y1)
    error = int(deltax / 2)
    y = y1
    ystep = None
    if y1 < y2:
        ystep = 1
    else:
        ystep = -1
    for x in range(x1, x2 + 1):
        if issteep:
            points.append((y, x))
        else:
            points.append((x, y))
        error -= deltay
        if error < 0:
            y += ystep
            error += deltax
    # Reverse the list if the coordinates were reversed
    if rev:
        points.reverse()
    return points

def computeHistogramEqual(buf):
    '''
    http://stackoverflow.com/questions/5858902/histogram-equalization-of-matplotlib-color-tables
    '''
    '''
    Non-linear colormap...
    http://stackoverflow.com/questions/22521382/nonlinear-colormap-matplotlib
    '''
    divs = 1000
    imvals = np.sort(buf.flatten())
    lo = imvals[0]
    hi = imvals[-1]
    steps = (imvals[::len(imvals)/divs] - lo) / (hi - lo)
    num_steps = float(len(steps))
    interps = [(s, idx/num_steps, idx/num_steps) for idx, s in enumerate(steps)]
    interps.append((1, 1, 1))
    cdict = {'red' : interps,
             'green' : interps,
             'blue' : interps}
    histeq_cmap = matplotlib.colors.LinearSegmentedColormap('HistEq', cdict)
    
    return histeq_cmap

class nlcmap(object):
    '''
    http://stackoverflow.com/questions/22521382/nonlinear-colormap-matplotlib
    '''
    def __init__(self, cmap, levels):
        self.cmap = cmap
        self.N = cmap.N
        self.monochrome = self.cmap.monochrome
        self.levels = np.asarray(levels, dtype='float64')
        self._x = self.levels
        self.levmax = self.levels.max()
        self.transformed_levels = np.linspace(0.0, self.levmax,
             len(self.levels))

    def __call__(self, xi, alpha=1.0, **kw):
        yi = np.interp(xi, self._x, self.transformed_levels)
        return self.cmap(yi / self.levmax, alpha)
    
def renderMultipleSlineSegments(toplist=None):
    
    '''
    plot line segments probability (additive) into buffer using
    http://en.wikipedia.org/wiki/Bresenham%27s_line_algorithm
    http://www.roguebasin.com/index.php?title=Bresenham%27s_Line_Algorithm
    '''
    
    buf = np.zeros( shape=im.shape )
    
    all_points = list()
    
    for slkey in sl_dict.keys():    
        sl = sl_dict[slkey]
        
        if toplist is not None and slkey in toplist:
            steps = len( sl[0] ) - 1
            for idx in range(0, steps):   
                #normalize by shifting / scaling to [0,1]
                prob = slines_prob[slkey][1][idx]    
                
                x0 = sl[0][idx]; x1 = sl[0][idx+1]
                y0 = sl[1][idx]; y1 = sl[1][idx+1]
                
                # scale to input LIC image buffer size
                x0sc = IMGSCALEX * x0
                y0sc = im.shape[1] - IMGSCALEY*y0
                x1sc = IMGSCALEX * x1
                y1sc = im.shape[1] - IMGSCALEY*y1
                
        #            points = get_line(int(x0),int(y0),int(x1),int(y1))
                points = get_line(int(round(x0sc)),\
                                  int(round(y0sc)),\
                                  int(round(x1sc)),\
                                  int(round(y1sc)))
                
                # add prob vel density to buffer
                for pt in points:
                    all_points.append(pt)
                    curr_prob = buf[pt[1], pt[0],0]
                    if prob > curr_prob:
                        buf[pt[1], pt[0],0] += prob
                    #print buf[pt[1], pt[0],0]
                
   
    #http://stackoverflow.com/questions/17201172/a-logarithmic-colorbar-in-matplotlib-scatter-plot
    '''
    The way that matplotlib does color mapping is in two steps, first a Normalize function 
    (wrapped up by the sub-classes of matplotlib.colors.Normalize) which maps the data you 
    hand in to [0, 1]. The second step maps values in [0,1] -> RGBA space.
    You just need to use the LogNorm normalization class, passed in with the norm kwarg.
    '''
    '''
    map = matplotlib.cm.ScalarMappable(norm=matplotlib.colors.LogNorm())
    map.set_array( buf[:,:,0] )
    map.set_cmap('gist_rainbow')#sequential
    map.autoscale()
    min, max = map.get_clim()
    '''
    min = np.amin(buf[:,:,0])
    max = np.amax(buf[:,:,0])
    
#    hmap = computeHistogramEqual(buf)

    rg = max-min
    
    #the choice of the levels depends on the data:
    levels = np.concatenate((
        np.linspace(max - 0.1*rg, max, 3),np.linspace(min, min + 0.1*rg, 3)
        ))
    
    levels = levels[levels <= max]
    levels.sort()
    
    cmap_nonlin = nlcmap(plt.cm.jet, levels)
    
    #for the colorbar we map the original colormap, not the nonlinear one:
    sm = plt.cm.ScalarMappable(cmap=plt.cm.jet, 
                    norm=plt.Normalize(vmin=min, vmax=max))
    sm._A = []
    
#    cbar = fig.colorbar(sm, cax=cbar_ax)
    #here we are relabel the linear colorbar ticks to match the nonlinear ticks
#    cbar.set_ticks(cmap_nonlin.transformed_levels)
#    cbar.set_ticklabels(["%.2f" % lev for lev in levels])
    
    fig, ax = plt.subplots()
    cbar = fig.colorbar(sm, orientation='horizontal')#, ticks=[min, max])
    cbar.set_ticks(cmap_nonlin.transformed_levels)
    cbar.set_ticklabels(["%.4f" % lev for lev in levels])
    plt.savefig('../img/colorbarCumulativeTest.png')
    
    #fig, ax = plt.subplots()
    #fig.colorbar(hmap, orientation='horizontal')#, ticks=[min, max] )
    #plt.savefig('../img/colorbarHistoCmap.png')
    
    maxp = np.amax(buf[:,:,0])
#    minp = np.amin(buf)
    
    for pt in all_points:
        #print x
        #print y
        #if buf[x,y,0] < 0.0001:
        #    buf[x,y,:] = im[x,y,:]
        #else:
        x = pt[1]; y = pt[0]
        #http://graph-plotter.cours-de-math.eu/
        #plotting this response for alpha
        #alpha = math.log( (buf[x,y,0] / maxp) + 1.0 ) + 0.3 
#        buf[x,y,:] = map.to_rgba(buf[x,y,0])
        buf[x,y,:] = cmap_nonlin(buf[x,y,0])       
#        buf[x,y,:] = map.to_rgba(hmap(buf[x,y,0])[0])
        #buf[x,y,3] = alpha

    # draw seed loc
    plt.figure() 
    plt.box('off')
    plt.imshow(im)
    imgplot = plt.imshow(buf) 
    
    plt.plot([ IMGSCALEX*sx ], [ im.shape[1] - IMGSCALEY*sy ], '+', color='black',linewidth = LW, ms=MS )
    
    plt.tick_params(\
    axis='both',          # changes apply to the x-axis
    which='both',      # both major and minor ticks are affected
    bottom='off',
    right='off',
    left='off',      # ticks along the bottom edge are off
    top='off',         # ticks along the top edge are off
    labelbottom='off', labelleft='off') # labels along the bottom edge are off
    
    plt.xlim([xmin, xmax])
    plt.ylim([ymin, ymax])
    
    plt.savefig('../img/cumulativeProbsSlinesTest.png')
    #plt.show()
    
def renderMultipleSlineSegments2(toplist=None):
    
    '''
    plot line segments probability (additive) into buffer using
    http://en.wikipedia.org/wiki/Bresenham%27s_line_algorithm
    http://www.roguebasin.com/index.php?title=Bresenham%27s_Line_Algorithm
    '''
    
    buf = np.zeros( shape=im.shape )
    
    all_points = list()
    
    for slkey in sl_dict.keys():    
        sl = sl_dict[slkey]
        
        if toplist is not None and slkey in toplist:
            steps = len( sl[0] ) - 1
            for idx in range(0, steps):   
                #normalize by shifting / scaling to [0,1]
                prob = slines_prob[slkey][1][idx]    
                
                x0 = sl[0][idx]; x1 = sl[0][idx+1]
                y0 = sl[1][idx]; y1 = sl[1][idx+1]
                
                # scale to input LIC image buffer size
                x0sc = IMGSCALEX * x0
                y0sc = im.shape[1] - IMGSCALEY*y0
                x1sc = IMGSCALEX * x1
                y1sc = im.shape[1] - IMGSCALEY*y1
                
                rr, cc, val = line_aa(int(round(y0sc)),\
                                  int(round(x0sc)),\
                                  int(round(y1sc)),\
                                  int(round(x1sc)))
                
                for idx in range(0, len(val)):
                    all_points.append((rr[idx], cc[idx]))
                    if buf[rr[idx], cc[idx], 0] < prob * val[idx]:
                        buf[rr[idx], cc[idx], :] = [prob * val[idx]]*4
                    
    min = np.amin(buf[:,:,0])
    max = np.amax(buf[:,:,0])

    rg = max-min
    
    #the choice of the levels depends on the data:
    levels = np.concatenate((
        np.linspace(max - 0.1*rg, max, 3),np.linspace(min, min + 0.1*rg, 3)
        ))
    
    levels = levels[levels <= max]
    levels.sort()
    
    cmap_nonlin = nlcmap(plt.cm.jet, levels)
    
    #for the colorbar we map the original colormap, not the nonlinear one:
    sm = plt.cm.ScalarMappable(cmap=plt.cm.jet, 
                    norm=plt.Normalize(vmin=min, vmax=max))
    sm._A = []
    
    fig, ax = plt.subplots()
    cbar = fig.colorbar(sm, orientation='horizontal')#, ticks=[min, max])
    cbar.set_ticks(cmap_nonlin.transformed_levels)
    cbar.set_ticklabels(["%.4f" % lev for lev in levels])
    plt.savefig('../img/colorbarCumulativeTestAA.png')
   
    maxp = np.amax(buf[:,:,0])

    for pt in all_points:
        x = pt[0]; y = pt[1]
        buf[x,y,:] = cmap_nonlin(buf[x,y,0])       
                
    # draw seed loc
    plt.figure() 
    plt.box('off')
    plt.imshow(im)
    imgplot = plt.imshow(buf) 
    
    plt.plot([ IMGSCALEX*sx ], [ im.shape[1] - IMGSCALEY*sy ], '+', color='black',linewidth = LW, ms=MS )
    
    plt.tick_params(\
    axis='both',          # changes apply to the x-axis
    which='both',      # both major and minor ticks are affected
    bottom='off',
    right='off',
    left='off',      # ticks along the bottom edge are off
    top='off',         # ticks along the top edge are off
    labelbottom='off', labelleft='off') # labels along the bottom edge are off
    
    plt.xlim([xmin, xmax])
    plt.ylim([ymin, ymax])
    
    plt.savefig('../img/cumulativeProbsSlinesTestAA.png')
    #plt.show()
    
def renderMultipleSlineSegments3(top=None):
    from matplotlib.collections import LineCollection
    
    toplist_dict = dict(toplist)
    
    fig = plt.figure()
    #ax = fig.gca(projection='3d')
    axes = fig.add_subplot(1, 1, 1, axisbg='white')
    #axes.set_aspect('equal', 'datalim')
    # render streamlines
    
    plt.box('off')
    plt.imshow(im)
   
    for slkey in sl_dict.keys():
        if top is not None and slkey in top:
            sl = sl_dict[slkey]
           
            try:
                l = list(slines_prob[slkey][1])
                #print 'length of list l: ' + str( len(l) )
                
                min_= min(l)
                max_ = max(l)
            except:
                #print 'streamline not more than one integration step'
                continue
            
            # draw each segment a color depending on prob density
            steps = len( sl[0] ) - 1
            
            x = list()
            y = list()
            lws = list()
            clrs = list()
            
            for idx in range(0, steps):
               
                #normalize by shifting / scaling to [0,1]
                scale = (slines_prob[slkey][1][idx] - min_0) / math.fabs( max_0 - min_0) 
            
                    
                prob = slines_prob[slkey][1][idx]
                print prob
                '''
                if scale < 0:
                    scale = 0.0
                print scale
                '''
                clr = list(mapper.to_rgba(x=prob)) 
                tot = toplist_dict[slkey]
                alpha = tot / max2 #normalize magnitude for alpha  
                clr[3] = alpha
                print clr  
                #clr = (1.0*scale, 1.0*scale, 1.0*scale, 1.0)
                x0 = sl[0][idx]; x1 = sl[0][idx+1]
                y0 = sl[1][idx]; y1 = sl[1][idx+1]
                
                # scale to input LIC image buffer size
                x0sc = IMGSCALEX * x0
                y0sc = im.shape[1] - IMGSCALEY*y0
                x1sc = IMGSCALEX * x1
                y1sc = im.shape[1] - IMGSCALEY*y1
                
               
                
                SCALE_W = 3.0
                MIN_W = 1.0
                x.append(x0sc); x.append(x1sc)
                y.append(y0sc), y.append(y1sc)
                lws.append(scale*SCALE_W + MIN_W)
                
                clrs.append(tuple(clr))
                
        points = np.array([x,y]).T.reshape(-1,1,2)
        segments = np.concatenate([points[:-1], points[1:]], axis=1)
        lc = LineCollection(segments, linewidths=lws, colors=clrs)
        axes.add_collection(lc)
#        plt.plot([ sx ], [ sy ], '+', color='black', linewidth=LW , ms=MS) # draw seed loc
    
    plt.plot([ IMGSCALEX*sx ], [ im.shape[1] - IMGSCALEY*sy ], '+', color='black',linewidth = LW, ms=MS )
    
    plt.tick_params(\
    axis='both',          # changes apply to the x-axis
    which='both',      # both major and minor ticks are affected
    bottom='off',
    right='off',
    left='off',      # ticks along the bottom edge are off
    top='off',         # ticks along the top edge are off
    labelbottom='off', labelleft='off') # labels along the bottom edge are off
    
    plt.xlim([xmin, xmax])
    plt.ylim([ymin, ymax])
    
    plt.savefig('../img/cumulativeProbsSlinesTest3.png')
    #plt.show()
    
if __name__ == '__main__':
    global prob, slines, MEM
    
    PERCENT = 0.01
    MEM = 990
    MAX_PER_FILE = 110
    sx = 60 # seed x coord
    sy = 60 # seed y coord
    
    #load background LIC
    im = plt.imread('/home/behollis/Dropbox/eurovis2015shortpaperCode/lockflowLICMeanField.png')
    IMGSCALEX = float(im.shape[0]) / 125. 
    IMGSCALEY = float(im.shape[1]) / 125.
    
    # register plots
    XSETm = 10
    XSETp = 40
    YSETm = 6
    YSETp = 28
    
    xmin = IMGSCALEX*sx - IMGSCALEX*XSETm
    xmax = IMGSCALEX*sx + IMGSCALEX*XSETp
    ymin = im.shape[1] - (IMGSCALEY*sy - IMGSCALEY*YSETm)
    ymax = im.shape[1] - (IMGSCALEY*sy + IMGSCALEY*YSETp)
     
    flist = openH5pyFiles()
    sl_dict = readStreamlines(sx, sy, b_files = flist, f_files = None, max_per_file=MAX_PER_FILE)
    closeH5pyFiles(flist)
    
    slines_prob, prob_list = calcStreamlineProbs2([sx, sy], sl_dict, MEM)
    
    mapper = matplotlib.cm.ScalarMappable()
    mapper2 = matplotlib.cm.ScalarMappable() #mapper for total prob streamlines
    min_0 = min(prob_list)
    max_0 = max(prob_list)
        
    shifted_prob = np.array(prob_list) #- min_0
    
    mapper.set_array( shifted_prob )
    mapper.set_cmap('jet')#sequential
    mapper.autoscale()
    
    #seed point
    LW = 10.0
    MS = 10.0
    
    # sort slines based on overall likelihood
    toplist = list()
    total_probs = list()
    for slk in sl_dict.keys():
        print slk
        toplist.append( (slk,slines_prob[slk][0]) )
        
    for key in slines_prob.keys(): 
        total_probs.append( slines_prob[key][0] )
    
    sorted( toplist, key=itemgetter(1) )
    
    toplist2 = toplist[0:int(len(toplist)*PERCENT)]
    topslines = list()
    for item in toplist2:
        topslines.append(item[0])
        
    #put most probable on top by reversing draw order with most prob last
    toplist2.reverse()
    topslines.reverse()
    
    mapper2 = matplotlib.cm.ScalarMappable() #mapper for total prob streamlines
    min2 = min(total_probs)
    max2 = max(total_probs)
    mapper2.set_array( np.array(total_probs) )
    mapper2.set_cmap('OrRd')#uential
    
    #norm = matplotlib.colors.Normalize(vmin=min_0, vmax=max_0)
    norm2 = matplotlib.colors.Normalize(vmin=min2, vmax=max2)
    mapper2.set_norm(norm2)
    #mapper.set_norm(norm)
    #mapper.autoscale()
    mapper2.autoscale()
    
    
    renderStandardPlot(topslines)
    renderColorbars()
#    renderOverallLikelihoods(topslines)
#    renderIndividualSlines()

#    renderMultipleSlineSegments(topslines)
#    renderMultipleSlineSegments2(topslines)
#    renderMultipleSlineSegments3(topslines)
   
    
    print 'finished!'
        